#include <stdio.h>
#include <stdlib.h>
#include <cv.h>
#include <highgui.h>
#include <math.h>


void CoocurrenceMatrix(IplImage *img, double **Cm, int nColor, int dX, int dY, int X, int Y) {
 int i,j;
 IplImage *img_quant;

    img_quant = cvCloneImage(img);
 
 // alocar uma nova imagem de tamanho maior para ser processada
 IplImage *img_b = cvCreateImage(cvSize(img_quant->width+((dY+1)*2), img_quant->height+((dX+1)*2) ), img_quant->depth, img_quant->nChannels);
 CvPoint offset = cvPoint(dX+1, dY+1);
 // copia a imagem quantizada e replica os pixels das bordas
 cvCopyMakeBorder(img_quant, img_b, offset, IPL_BORDER_REPLICATE, cvScalarAll(0));
 
 uchar *data = (uchar *)img_b->imageData;
 int step = img_b->widthStep;
 int height = img_b->height;
 int width = img_b->width;
 
if (X == 0){
 // percorre a imagem
 for (i = dX; i < height-dX; i++) {
      for (j = dY; j < width-dY; j++) {
	  int pref = data[i*step+j]; //pixel referencia
	  int pviz = data[(i+dX)*step+(j+dY)]; //pixel vizinho
	  Cm[pref][pviz]++;
      }
 }
}
else {
 // percorre a imagem
 for (i = X; i < X + 128 -dX; i++) {
      for (j = Y; j < Y + 128 -dY; j++) {
	  int pref = data[i*step+j]; //pixel referencia
	  int pviz = data[(i+dX)*step+(j+dY)]; //pixel vizinho
	  Cm[pref][pviz]++;
      }
 }
}
 double sum = 0;
 // calcula a soma das ocorrencias
 for (i=0; i<nColor; i++) {
      for (j=0; j<nColor; j++) {
	    sum+= Cm[i][j];
      }
 }
 // normaliza a matriz de forma que sum(Cm) = 1
 for (i=0; i<nColor; i++) {
      for (j=0; j<nColor; j++) {
	    Cm[i][j] /= sum;
      }
 }
}

void Haralick(double **Cm, int size, double *HDesc) {
    
    int i,j;
    double m_r = 0.0;
    double m_c = 0.0;
    double s_r = 0.0;
    double s_c = 0.0;
    double *Pi = (double *) calloc(size, sizeof(double));
    double *Pj = (double *) calloc(size, sizeof(double));
  
    for (i = 0; i<size; i++) {
	for (j = 0; j<size; j++) {
	      Pi[i] += Cm[i][j];
	}
	m_r += i*Pi[i];	
    }
    for (j = 0; j<size; j++) {
	for (i = 0; i<size; i++) {
	      Pj[j] += Cm[i][j];
	}
	m_c += j*Pj[j];	
    }
    
    for (i = 0; i<size; i++) {
	s_r += ((i-m_r)*(i-m_r)) * Pi[i];
	s_c += ((i-m_c)*(i-m_c)) * Pj[i];
    }
    s_r = sqrt(s_r);
    s_c = sqrt(s_c);
    
    double maxp = 0.0;
    double corr = 0.0;
    double cont = 0.0;
    double unif = 0.0;
    double homo = 0.0;
    double entr = 0.0;
    double auxv = 0.0;
    
    for (i = 0; i<size; i++) {
      for (j = 0; j<size; j++) {
	  auxv = Cm[i][j];
	  // Descritor 1 - Maxima Probabilidade
	  if (maxp < auxv) {
	      maxp = auxv;
	  }
	  // Descritor 2 - correlacao
	  if (s_r > 0 && s_c > 0) {
	      corr+= ((i-m_r)*(j-m_c)*auxv) / (s_r*s_c);	    
	  }
	  // Descritor 3 - contraste
	  cont += ( (i-j)*(i-j)*auxv );
	  // Descritor 4 - energia (uniformidade)
	  unif += (auxv*auxv);
	  // Descritor 5 - homogeneidade 
	  unif += (auxv) / (1 + abs(i-j));
	  // Descritor 6 - entropia
	  if (auxv != 0) {
	    //entr+= auxv*( log(auxv) / log(2) );
	    entr+= auxv*( log2(auxv));
	  }
      }
    }
	
    entr = -entr;
    
    HDesc[0] = maxp;
    HDesc[1] = corr;
    HDesc[2] = cont;
    HDesc[3] = unif;
    HDesc[4] = homo;
    HDesc[5] = entr;  
}

double distManhattan(double *p, double *q, int size) {
    int i;
    double dist =0;
    for (i=0; i<size; i++) {
	  dist += fabs(p[i]-q[i]);      
    }
    return dist;  
}
